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ABSTRACT 

We compare the structure of star-forming molecular clouds in different regions of Orion A to determine how the column density 
probability distribution function (A-PDF) varies with environmental conditions such as the fraction of young protostars. A corre¬ 
lation between the A-PDF slope and Class 0 protostar fraction has been previously observed in a low-mass star-formation region 
(Perseus) by Sadavoy; here we test if a similar correlation is observed in a high-mass star-forming region. We use Herschel PACS 
and SPIRE cold dust emission observations to derive a column density map of Orion A. We use the Herschel Orion Protostar Survey 
(HOPS) catalog for accurate identification and classification of the Orion A young stellar object (YSO) content, including the cold and 
relatively short-lived Class 0 protostars (with a ''0.14 Myr lifetime). We divide Orion A into eight independent 0.25 square degree 
(13.5 pc^) regions; in each region we fit the A-PDF distribution with a power-law, and we measure the fraction of Class 0 protostars. 
We use a maximum likelihood method to measure the A-PDF power-law index without binning the column density data. We find that 
the Class 0 fraction is higher in regions with flatter column density distributions. We test the effects of incompleteness, extinction- 
driven misclassification of Class 0 sources, resolution, and adopted pixel-scales. We show that these effects cannot account for the 
observed trend. Our observations demonstrate an association between the slope of the power-law A-PDF and the Class 0 fractions 
within Orion A. Various interpretations are discussed including timescales based on the Class 0 protostar fraction assuming a constant 
star-formation rate. The observed relation suggests that the A-PDF can be related to an “evolutionary state” of the gas. If universal, 
such a relation permits an evaluation of the evolutionary state from the A-PDF power-law index at much greater distances than those 
accesible with protostar counts. 
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1. Introduction 


2. Column density maps and source catalogs 


The structure of molecular clouds can be described by the prob¬ 
ability di stribution function (A- PDF) of their colu mn densi¬ 
ties (e.g., Kainula inen et alJl2009l) . Observations (e .g.. Hill et al 


2011 


[Hughes et alJ 2013 ) and theoretical studies (iPadoan et al 


20131) suggest that the A-PDF dep ends on environmen tal con¬ 
ditions s uch as turbulence ( e.g., iPadoan et al.l Il997h . grav- 
ity (e.g., IGessen^et al.i l2000t) . magnetic field strength (e.g., 
Molina et al.ll201 2[). and star-formation (e.g., Schneider et al.l 
2013l:lKainulainen et al.ll20T4t lAbreu et alJl2015l) . In particular. 


an increasing A-PDF slope has been shown to be correlated with 
an increasing fr action of Class 0 protostars in individual clumps 
within Perseus (ISadavovll2013h . 


Here we aim to test if a similar correlation is observed in a 
different star-forming region. We present a study of subregions 
within the Orion A molecular cloud that quantifies the link be¬ 
tween the incidence of Class 0 protostars and the A-PDF slope. 
We bring together two data sets: column density measurements 
using Herschel observations, and an accurate protost ellar census 
from the Herschel Orion Protostar Survey (HOPS; IStutz et ^ 
[20131: [Fischer et al.[[2013l Furlan et al., in prep). With this anal¬ 
ysis, we quantify spatial variations in the relationship between 
column density distributions and fractions of young protostars 
within Orion A. 


Send offprint requests to: stutz@mpia.de 


2.1. Column density maps 

We present the Orion A column density, N(H), map (Figs. [1] 
and lA.ll) derived from Herschel 160 pm to 500 pm emission 
maps calibrated against Planck and IRAS data ([Bernard et al.[ 
[20101) . The da ta were observed as pa rt of the Herschel Gould 
Belt program ([Polvchroni et al.l[2013[) . The column dens i ty and 
tem perature maps were de rived as in [Stutz et ^ ([20101. [2013[) 
and lLaunhardt et akl ([20131) . See Appendix [Al for more details. 


2.2. Catalog of Class 0, Class I, and flat spectrum young 
stellar objects 

The young stellar object (YSO) catalog is the union of 
the PACS Bright Red Sources (PBRS) sample (IStutz et al.[ 
[2013[) of extremely young Herschel-detected Class 0 p roto- 
stars ([Heiderman & Evansll2015[ : [van Kempen et al.l[20091) and 
the HOPS sample (Furlan et al., in prep.. and iFischer et alJ2013[) 
of Class 0, Class I, and flat spectrum YSO s. The HOPS YSO cat - 
alog is based on the Spitzer catalog from [Megeath et aP ([20121) . 
which excludes extragalactic and stellar source contamination 
by means of infrared colors. The YSOs have well-sampled 
SEDs from near-infrared to the submillimeter wavelengths, in¬ 
cluding Herschel 70 pm, 100 pm, and 160 pm measurements 
as well as APEX LABOCA 870 pm data. We also include 
APEX SABOCA 350 pm data when available ([Stanke et al.[ 
[ 2 OIOI : [Stutz et al][2013t [Safron et^[2015l) to sample the peak 
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Fig. 1. Orion A N(H) map shown on a log scale; the N(H)= 
3.7 timeslQ^^ cm~^ contour is indicated in grey. The locations 
of protostars are indicated as x-symbols. Numbered boxes (with 
1/4 square degree area or 3.67 pc on a side) indicate the regions 
into which we divide the Orion A cloud. 

of the cold envelope emission. The locations of the YSOs are 
indicated in Fig. [1] 

3. Results 

3.1. Power-law N(H) distributions 

We analyze the high column density portion of the A-PDFs 
within the filamentary sub-regions of Orion A. iLombardi et al.l 
(l20L5ll demonstrate that the high-N(H) portion of the A-PDF 
can be reliably inferred from Herschel observations above Ak ~ 
0.2 mag even when accounting for complicating factors such as 
line-of-sight blending and errors in the background emission or 
extinction estimates. We divide Orion A into 8 separate square 
regions that are 3.67 pc on a side. The region locations are cho¬ 
sen such that they are independent (have no overlap), avoid the 
center of the Orion Nebula Cluster (ONC) which is affected by 
gross protostellar incompleteness and saturation in both Spitzer 
and Herschel observations), and cover regions that are presently 
forming stars. The region locations are shown in Fig.[T] 

We extract the A-PDF and ht the high column density por¬ 
tion with a power-law in each region. Our results are presented 
in Table [1] and Figs. |2] and [3] The Herschel 250 pm beam is 
~ 18" FWHM; therefore we use an 18"pixel scale for extracting 
the N(H) distributions to minimize pixel-to-pixel correlations. 
Inspection of the A-PDFs reveals that they have an approximate 



log N(H) [cm"^] 

Fig. 2. N(FI) distributions for each region shown in Fig. [1] The 
combined N(H) distribution for all regions is shown in grey. The 
dashed curves show the best-fit slope assuming a power-law dis¬ 
tribution of N(H); each value for the slope is listed (see Table [1] 
for slope errors). Here we show PDF data of bins containing 
more than 10 pixels, and each major tick mark represents a fac¬ 
tor of 10. Mean-normalization, N(H)/<N(H)> in each region, 
has no effect on the derived PDF slopes because the operation 
simply results in a linear translation of the x-axis. 


power-law shape above log N(H) ~ 22.0 to 22.4 and below log 
N(H) ~ 22.7 to 23.0. Thus, for each region the definition of the 
fitting regime is driven by the requirement of excluding areas 
with curvature in the A-PDFs. The minimum and maximum 
N(H) values used to derive the power-law indices are listed in 
Tabled 

We fit a power-law to the A-PDF using a maximum like¬ 
lihood method that does not require binning the data and is 
therefore free of any potential errors due to bin size. We obtain 
power-law indexes that vary from a - -0.93 to a = -2.95, or a 
factor of about 3. We test the effects of varying the beam size and 
pixel size on our derived indices and find that our fitting method 
is robust. We also test the effect of varying the upper bound of 
the fit and find that the effect on the indices is much smaller than 
the measured variations across regions. See appendixlBlfor more 
details. 
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Table 1. Best-fit A-PDF power-law slopes and Class 0 protostar fractions 


Region" 

R.A. 

Decl. 

[Min,Max] log NIH)* 

PDF index" 

NciassO 

Nyso 

Fraction" 

1 

05:35:11.91 

-05:01:08.00 

[22.0, 22.9] 

-0.93 ± 0.04 

24 

45 

0.53 ± 0.07 

2 

05:35:07.37 

-05:43:23.15 

[22.0, 22.8] 

-1.03 ±0.04 

13 

27 

0.48 ±0.10 

3 

05:35:44.43 

-06:13:49.45 

[22.2, 23.0] 

-1.95 ±0.05 

8 

22 

0.36 ±0.10 

4 

05:36:36.13 

-06:43:45.72 

[22.1,22.9] 

-2.26 ± 0.06 

9 

16 

0.56 ±0.12 

5 

05:38:38.41 

-07:01:35.38 

[22.1,22.7] 

-2.04 ± 0.06 

7 

30 

0.23 ± 0.08 

6 

05:39:46.66 

-07:31:29.36 

[22.4, 23.0] 

-2.30 ± 0.09 

7 

26 

0.27 ± 0.09 

7 

05:41:07.68 

-08:02:05.33 

[22.2, 23.0] 

-1.91 ±0.05 

8 

48 

0.17 ±0.05 

8 

05:42:23.95 

-08:41:22.42 

[22.1,22.7] 

-2.95 ± 0.08 

3 

20 

0.15 ±0.08 

AIF 



[22.2, 23.0] 

-1.90 ±0.02 

79 

234 

0.34 ± 0.03 


“ Region numbers as in Fig.[TJ each region is 0.5 deg in a side. 

Minimum and maximum log N(FI) values respectively included in the fit to the power-law index of the PDF. 

PDF power-law index errors are estimated using the = 1.0 interval. 

Fraction of Class 0 protostars relative to the total number YSOs: Nciasso/Nyso. The fractional errors are derived assuming a binomial distribution 
(see text). 

'■ Integrated properties of the above regions. 


3.2. Class 0 protostar fractions 

Using the catalogs discussed in § 12.21 we count the total num¬ 
ber of YSOs and the subset of Class 0 protostars in each re¬ 
gion (see Table [1] and Fig. [T]i. We define the Class 0 fraction 
as the number Class 0 protostars divided by the total number 
of YSOs: NciassO / Nyso- We dehne the Class 0 protostars as 
the subset of sour ces with bolometric temperatures Tboi < 70 K 
(IChen et al.lll995h . Tboi is dehned as the temperature of a black- 
body with t he same flux-weighte d mean frequency as the ob¬ 
served SED dMvers & Laddll99^ . The original protostellar Tboi 
classihcation acco unted for the effects of foreground extinction 
dChen et al ; see also Dunham et al. 2013). Here we calcu¬ 

late Tboi based on the observed SED, without additional extinc¬ 
tion corrections. See text below and AppendixIClfor analysis on 
the effects of foreground extinction. Our wavelength coverage 
samples the peak of even the coldest Class 0 protostellar SEDs 
and allows for a robust protostellar cl assihcation based on Tboi 
dStutz et al.ll2013l : iDunham et al.ll2014b . 

We estimate the errors on the fraction of Class 0 protostars 
to the total number of YSOs using the binomial distribution as 
follows. In each region the number of Class I protostars is (ob¬ 
viously) Nciassi = Nyso“ NciassO, where we include the flat- 
spectrum YSOs in the Class I sample. Since the YSOs in each 
region have a probability p of being Class I and a probability q 
= (1 - p) of being Class 0, the expected number of Class 0 proto¬ 
stars is therefore Nyso R Q- The NciassO error is V-^yso x px q 
and the fractional error is therefore cr = The hnal 

numbers of YSOs, Class 0 protostars, fractions and respective 
errors in each region are presented in Table [1] 

Two principal effects could potentially alter the measured 
fractions of Class 0 protostars: incompleteness and misclassih- 
cation. Variations in flux completeness across Orion A are dom¬ 
inated by the spatially-varying level of nebulosity. We deter¬ 
mine the completeness limits across the cloud by injecting fake 
sources into the PACS 70 pm images and measuring the flux at 
which 90% of the sources are recovered. Region 1 has the high¬ 
est 70 pm flux completeness limit. We apply this limit to the 
entire YSO sample in each region and And that the Class 0 frac¬ 
tions are largely unaffected by incompleteness. We therefore do 
not correct the observational numbers presented in Table [T] See 
AppendixICland Eig. IC.ll for more discussion. 

With a reliable sample of YSOs and protostars, the main 
classification ambiguities that may hinder the identification of 



N(H) slope a 


Fig. 3. Correlation between the power-law index of the A-PDE 
and Class 0 fraction across Orion A regions. The Class 0 frac¬ 
tion can be related to a timescale as indicated assuming a con¬ 
stant SER (see text). Grey x-symbols are th e data from Perseus 
from ISadavovI (l2013h , ISadavov et al.l (l2014l) . and Sadavoy (pri¬ 
vate communication, 2015). Dashed grey curve shows the linear 
fit to the data, which have a correlation coefficient of 0.7. 


Class 0 protostars include any extrinsic effects (that is, effects 
acting independently of the Class 0 envelope) that will cause the 
observed SED to appear redder. The two main cuprits are the 
inclination of the disk relative to the line-of-sight (LOS) and 
foreground reddening. Using a grid of pr otostellar SED models 
(Furlan et al., in prep., IStutz et ani2013h . we test the effects of 
both on the Tboi-hased classification. Even assuming the most 
elevated levels of extinction measured from our N(H) maps to¬ 
ward the positions of the protostars, we And that neither effect 
can account for the variations in the Class 0 fractions reported in 
§ 13.31 See AppendixICland Eigs. IC.3l and lC.4l for more details. 
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3.3. The connection between the N-PDF and the Class 0 
fraction 

The observed A-PDF slopes and Class 0 fractions show a spatial 
variation within Orion A. From North to South, the A-PDF slope 
steepens and the C lass 0 fra c tion d e creases. Fig. which also 
includes data from ISadavovI (l2013h . ISadavov et all(l2014h . and 
Sadavoy (private communication, 2015), shows this observed 
correlation. A linear fit to the data results in a slope of 0.2 and 
an intercept of 0.7. We measure a correlation coefficient of 0.7. 
According to the t-test there is a 0.3 % chance of having a cor¬ 
relation coefficient this high in an uncorrelated random sample. 

4. Discussion 

The observed A-PDF slopes and Class 0 fractions indicate clear 
spatial variations within a single molecular cloud (Orion A). 
Making the simple assumption of a constant star-formation rate 
(SFR), the Class 0 fraction can be related to time (t). If the 
Class 0 and Class 1 lifetimes are 0.1 Myr and 0.4 Myr, re¬ 
spectively, then at f < 0.1 Myr the fraction will be equal to 
1; at 0.1 Myr > f > 0.5 Myr the fraction is = [f/(0.1Myr)]“'; 
at f > 0.5 Myr, the fraction will be constant and equal to 
1/(1 -H fi/fo) ~0.26, where to and fi ar e the lifetimes of th e 
Class 0 and Class 1 phases, respectively (iDunham et al.l 120141) . 
Under these assumptions, high fractions correspond to short evo¬ 
lutionary timescales of the on-going star-formation event, and 
the A-PDFs change on ~ 0.3 Myr timescales. 

However, variations in the SFR could explain the observed 
trend. We kn ow that Orion A has formed stars for longer than 
0.5 Myr (e.g.. lMegeath et aD2012h . An elevated Class 0 fraction 
may therefore reflect an increasing SFR over time compared to 
regions with lower Class 0 protostar fractions. 

O-star feedback could compress the gas and cause a flat¬ 
tening of the A-PDF slopes. Regions 1 and 2 would be most af¬ 
fected by feedback because of the O-star population in the north¬ 
ern portion of Orion A. In Perseus O-stars themselves are not 
responsib le because the r egion does not contain any such stars. 
However, ISadavov et al.l (l2014h propose that feedback from the 
low-mass protostars themselves may account for the observed 
correlation between Class 0 fraction and A-PDF slope. This 
mechanism could potentially operate in Orion as well. 

5. Conclusions 

We analyze the A-PDF slopes and Class 0 protostellar fractions 
within sub-regions of the Orion A molecular cloud. Our conclu¬ 
sions are as follows. 

• We observe a progression from shallow to steep A-PDF slopes 
from North to South in the Orion A cloud. This progression 
shows that there is no unique A-PDF, but that the A-PDF shape 
depends on environment. 

• We observe a correlation between increasing A-PDF slope 
and increasing fraction of Class 0 protostars in subregions of 
Orion A. Under the assumption that the Class 0 fraction is re¬ 
lated to time via an assumption of a constant SFR, evolutionary 
timescales for each region can be derived. This suggests that re¬ 
gions with shallower slopes have younger “evolutionary states”. 

If universal, this relation permits an evaluation of the evo¬ 
lutionary state from the N(H) power-law index measurement, 
which is possible at much greater distances than regions that are 
accesible with protostar counts. A key aspect of this study is that 
the slopes do not change significantly with resolution; therefore 


for fixed angular resolution we expect to obtain the same slope 
measurement over a broad range of distances. 
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Appendix A: Column density map of Orion A 

We have generated reduced data products for the Orion A re¬ 
gion using HIPE processing to level 1 followed by final level 2 
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Scanamo rphos processi ng (version 24.0, using the “galactic” 
option, iRoussell l2013h . The column density (and tempera- 
ture) maps were derived in a similar way as those pr esented in 
IStutz et ^ (I 2 OIOI l2013h . and iLaunhardt et alJ (1201 3h . with the 
total power emission levels at each of the four wavelength s cal¬ 
culated using Planck and IRAS data (iBernard et al.ll2()l^ . We 
briefly summarize the steps we use to generate the N(H) and 
temperature maps here and refer the reader the above works for 
more details. We apply the total power emission level correc¬ 
tions for each of the 160 /rm, 250 um, 350 fim, and 500 /tm 
intensity maps (IBernard et alj|20fo|) . We convolve the data to 
the largest beam; the 500 beam, with ~ 38" F WHM. We use 
the azi muthally averaged convolution kernels from lAniano et alJ 
(I 2 OI Ih . We then re-grid the data to a matched coordinate sys¬ 
tem and pixel scale, in this case we adopt an 18"pixel scale (see 
Appendix iBl for more details). The conversion to matched units 
assumes the beam sizes listed in the SPIRE instrument hand¬ 
book. We fit the spectral energy distribution (SED) of each pixel 
assuming a modified black-body function of the form 

5v = QBv(v,rd)(l-e-"'’'^), (A.l) 

where Q. is the solid angle of the emitting element, By{TS) is 
the Planck function at a dust temperature Ta, and t(v) is the op¬ 
tical depth at frequency v. Here, the optical depth is given by 
t(v) = Ah Rg] k{v), where Ah — 2 x A(H 2 ) -i- A(H) is the 
total hydrogen column density, mn in the proton mass, Ky is the 
assumed frequency dependent dust o pacity law, and R^d i s the 
gas-to -dust ratio, assumed to be 110 dSodroski et al.l[T99^ . We 
use the lOssenkopf & Hennind (1 19941) model dust opacities cor¬ 
responding to column 5 of their Table 1 (sometimes referred to as 
the “OH5” opacities). These opacities are meant to reflect grains 
with thin ice mantles after 10^ years of coagulation time at an 
assumed gas density o f 10^ cm“^. See I Stutz et ^ (1201 3h and 
ILaunhardt et all (1201 3l) for discussions on the systematic uncer¬ 
tainties introduced by the model dust opacity assumption. In or¬ 
der to apply the color and beam size corrections recommended in 
the SPIRE and PACS instrument handbooks, we first fit the un¬ 
corrected pixel SED to estimate the temperature. We then apply 
the interpolated correction value at that temperature and re-fit 
the SED. 

In an additional step, we use the 500 jim resolution temper¬ 
ature iji) map to convert the 250 jum intensity map to a column 
density map in order to improve the final resolution (Eigs. lA.ll 
and[T]i. Both maps compare well; only minor differences are ap¬ 
parent on the smallest scales caused by resolution effects, as ex¬ 
pected. The previously published Herschel-dt nvtA N(H) maps 
compare well to the maps we present here ('e.g.. lLombardi et^ 
120141: iPolvchroni et al.ll2013l) . 

Appendix B: Fitting N(H) PDFs without binning 


Table B.l. Power-law index error comparison 


Region 

Index 

error 

(T^ 

1 

-0.93 

0.039 

0.042 

2 

-1.03 

0.036 

0.035 

3 

-1.95 

0.047 

0.038 

4 

-2.26 

0.062 

0.053 

5 

-2.04 

0.061 

0.044 

6 

-2.30 

0.090 

0.067 

7 

-1.91 

0.047 

0.038 

8 

-2.95 

0.080 

0.067 


The errors reported here correspond to those derived from the N(H) 
intervals listed in Table[T] “ Error from eauation lB.il 


It is straightforward to show that the likelihood is maximized 
when the models are restricted to those with Uexp = Wobs- Hence, 
we simplify notation by n = riexp = ^obs and obtain 

n 

In X = a ^ In Z, H- n(ln v - + In n) 

i~\ 

To find where X is maximized, we differentiate w.r.t. a (or x = 
a -H 1) and set to zero 


Z n 

lnZ, + -- 


^max I*' 


•Z^. lnZ„ 

mm “ 


Z x _ yx 

mnY ^ 


which can be rewritten as 


where 


(InZ) = InZn^ 


<lnZ> = iy lnZ„ 

n 4-H 


rW - 1 1 


R' — 1 X ’ 


In Zniax _ ^max 

K = -. 


InZn, 


The eiTor is given bv lGouldl(ll995l) Eq. 2.4; 
1 


o-(x) 


1 R-^(lnR)2 

y" (R-- -1)2 


which can be written more elegantly as 
1 1 


o-(x) |x| 


1 


Q 


sinh Q 


<2 = 


xlnR 


Here we describe the method we use for fitting an index cr to a 
power-law distribution of the form dNIdZ oc Z" over a finite 
range in Z values from Zmin to Z^ax. The differential probability 
of column density Z,, given a power-law distribution Z" and total 
number of expected detections in the interval between Z^in 
and Zmax is 

Pi - Zf- - -Hexp- {x = a +1) 

I I ryx _ yx ' 

•^max -^min 

The likelihood X is given by the product of the individual prob¬ 
abilities, or 

^obs ^obs 

InX = = a^lnZ,+«ob.s(ln;c-ln(Z^3^-Z^J+ln?7exp)- 

/=i 1=1 


This expression can be Taylor expanded; 



(xlnR)2 

40 



(B.l) 


Eauation IB. 1 [ represents the analytical error solution assum¬ 
ing the power-law model accurately reflects the distribution of 
data values. Therefore it represents a lower limit for the errors. 
Alternatively, the errors can be estimated using the = 1.0 
interval, which compares well to those derived using the mini¬ 
mum variance bound (see Table ED for a comparison between 
the two). The power-law indices and = 1.0 errors are listed 
in Table [U 
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Fig.A.l. Orion A N(H) map, shown on a log scale. 
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Fig. B.l. Top; Best-fit power-law index for each region as a func¬ 
tion of N(H) map pixel size. The slope values are not affected by 
correlated pixel, and differences are 2% at most. Bottom: Error 
for each region as a function of N(H) map pixel size. The error 
increases with pixel size because number of counts decreases. 
All quantities shown here have been derived using the same 
N(H) limits presented in Table [1] 


B. 1. The effects of resolution, pixel-size, and fitting range on 

the power-law index 

As described above, the power-law indices for each region are 
extracted from 18" pixel N(H) maps derived from the 250 pm 
map of Orion A. Here we test the effect of adopting different 
pixel scales for the 250 pm N(H) map. As shown in Fig. IB. li the 
adopted pixel scale has a negligible effect on the indices, with 
a maximum effect on the best-fit indices of ~ 2%. The frac¬ 
tional errors increase with pixel size because the number of pix¬ 
els decreases. We find similar results using the 500 pm N(H) map 
(with a beam size of ~ 36" FWHM) and adopting pixel sizes of 
10", 20", 30", and 40"; the power-law index is only marginally 
affected by the choice in pixel size. The fractional errors exhibit 
the same behavior as for the 250 pm, but are somewhat larger 
(maximum of 10% at 40") due to the smaller number of pixels 
available to fit. Finally, in Fig. IB.2l we compare the indices de¬ 
rived from the 250 pm 18" pixel map with those derive from 
the 500 pm 40" pixel map. We find good agreement between the 
two slope estimates. 

We also test varying the upper bound of the fitting range 
shown as Max log N(H) in Table[T] For a change +A log N(H) = 


Fig.B.2. Indices derived from the 250 pm 18"pixel map vs. 
those derived from the 500 pm 40"pixel map, both using the 
N(H) limits presented in Table[T] 


0.3 we find the indices change by ~ 10% or less, a variation that 
is much smaller than differences between regions. 


Appendix C: The effects of completeness and 
extinction on the Class 0 protostar fraction 

C. 1. Completeness of protostellar catalogs 

The northern portion of Orion A (regions 1 and 2) are those most 
affected by incompleteness due to elevated levels of nebulosity. 
We estimate the completeness limit in the HOPS protostellar 
catalog as follows. We inject artificial sources into the 70 pm 
images and recov er them with our source-finding software (see 
IStutz et ^120131 for details). We estimate the mean 90% com¬ 
pleteness level for regions 1 and 2 to be 0.12 Jy, while for the 
rest of F1641 it is 0.03 Jy. We therefore apply a 0.12 Jy 70 pm 
flux limit to the HOPS protostar catalog, eliminating about 10% 
of protostars from the catalog. We then calculate Class 0 frac¬ 
tions from the remaining sources. We find excellent agreement 
between the completeness corrected and total raw sample frac¬ 
tions, as shown in Fig. 1C.II 

C.2. The effect of foreground extinction on Tboi 

Our goal is to assess the effects of foreground extinction on the 
Tboi-hased YSO classification. We assume here that the main 
contribution to an erroneous Tboi classification is the misidenti- 
fication of Class I sources as Class 0 protostars (Tboi) < 70 K). 
While other sources of contamination are also possible, scrutiny 
of the HOPS protostellar SEDs reveals that the current classifi¬ 
cation based on observed colors and fluxes, in combination with 
PACS and 870 pm data, is robust (Furlan et al., in preparation). 

We obtain observational constraints on the levels of 
foreground extinction toward individual protostars from the 
Herschel N(H) map. We measure the N(H) toward each pro¬ 
tostar using an annulus of size 18"x[1.5,4] (or ~ 11400 AU 
to 30300 AU), adopted to avoid extinction intrinsic to the pro¬ 
tostellar envelopes. Within each annulus we calculate the me¬ 
dian value of N(H). In Table 1C. II we list the median, mini¬ 
mum, and maximum protostellar N(H) values for each region, 
assuming a conver sion of Av/N(H)- 1.85 x 10^'magcm^^ (e.g., 
iBohlin et al.lll978l , derived for the diffuse ISM, and used here 
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Fig.C.l. Class 0 fraction in regions. Diamonds indicate the to¬ 
tal fraction (Table [TJ, while squares indicate the completeness 
corrected values. Results remain consistent with and without the 
application of the completeness limit. 


Table C.l. Upper limits to the distribution of protostellar extinc¬ 
tion per region 


Region 

Median" 

Ay 

Minimum" 

Ay 

Maximum" 

Ay 

1 

23.3 

0.64 

35.6 

2 

13.3 

2.18 

30.0 

3 

16.8 

4.34 

31.8 

4 

7.51 

4.08 

26.4 

5 

13.0 

2.63 

21.2 

6 

14.9 

4.03 

29.7 

7 

12.5 

3.65 

35.1 

8 

8.79 

4.54 

16.2 


“ Median, minimum, and maximum Ay (mag) from the calculated from 
the Herschel N(H) maps centered on the positions of protostars in each 
region. Ay are calculated in an annulus of size 18"x[1.5, ^pixels] (or ~ 
11400 AU to 30300 AU). 



Fig. C.l. Model Tboi values versus SED extincted Tboi.ext values. 
Colors indicate the assumed levels of extinction. Here we show 
500 randomly chosen models. As the extinction Ay increases, 
the model Tboi.ext decreases and the levels of contamination at 
Tboi.ext < 70 K increases. 


for notational simplicity). These values are upper limits because 
the N(H) map integrates the extinction through the entire cloud. 
The protostellar Ay values have a maximum of 35 mag in re¬ 
gions 1 and 7. However, more typical values range from ~20 
mag to ~10 mag. 


In order to assess the levels of contamination in the Class 
0 sample we us e the protostellar model grid presented in 
IStutz et ^ d2013l) . We refer the reader to that publication and to 
Furlan et al. (in prep.) for more details. In brief, we vary 5 param¬ 
eters in our grid; inclination of the disk relative to the LOS, the 
envelope density, the cavity opening angle, the disk size, and the 
luminosity of the central protost ar. Our grid contain s uniformly 
sampled paramaters and uses the lOrmel et akl (1201 ll) model dust 
opacities (“icsgra3”). In total, our grid contains 30400 individual 
models. About ~30% of the models have Tboi <70 K, similar to 
the observed protostellar distribution. We attenuate each model 
SED (adopting the wavelength coverage of the observations) 
with a range of Ay values between 1 and 35 mag, spaning the 
observed range in Ay toward the protostar sample (Table ICTt . 
For consistency with the above protostellar YSO measurements, 
we attenuate the model SEDs with the same model a s that used 
to derive the N(H) map (lOssenkopf & Hennindll99^ ; “OH5”), 
assuming Ay/AK = 14 (iNielbock et ak 20 12^ . In Fig. 1C.21 we 
show resulting extincted Tboi.ext values as a function of Tboi for a 
random subset of models. We then measure the fraction of mod¬ 
els that have Tboi > 70 K and Tboi.ext < 70 K as a function of Tboi 
and Ay, marginalizing over a uniform Ay distribution. The prob¬ 
ability that a given Class I protostar is masquerading as Class 0 
protostar is shown in Fig. 1C.31 top panel. The bottom panel of 
Fig. 1C.31 shows the probability that a given Class 0 protostar is 
correctly identified (has Tboi <70 K). 


We obtain a maximum fraction of ~ 9% of Class I sources 
that may appear as Class 0 protostars if models are extincted 
by up to Ay = 35 mag. For all models with Tboi >400 K, there 
is zero contamination of Class 0 sources up to extinctions of 
Ay ~ 60 mag. We conclude that the main source of extinction- 
driven contamination on the Class 0 sample is therefore from 
Class I sources with 70 K<Tboi <400 K. We also find that there 
is < 5% contamination of Class 0 protostars below Tboi.ext ~45 K 
and below foreground extinction levels of Ay = 35 mag. For a 
median extinction value of Ay = 20 mag, the fraction of Class 0 
contamination is less than 5% for Tboi.ext <55 K. 

Using the total contamination fractions shown in Fig. 1C.31 
and both the median and maximum observed Ay levels 
(Table ICTTi. we apply correction factors to the observed num¬ 
bers of Class 0 protostars, keeping the total number of YSOs 
fixed at the observed value. The dependence of the Class 0 frac¬ 
tion on extinction contamination is presented in Fig. 1C.41 While 
the overall amplitude of the fraction is affected by the foreground 
extinction effects, the shape is not. As noted above, the bottom 
panel of Fig. |C3] demonstrates that at Tboi.ext <55 K the amount 
of Class 0 contamination is less than 5% for our expected levels 
of extinction. We therefore test the effects of adopting a 55 K 
division between Class 0 and Class I protostars. We find that 
the results are similar to those shown in Fig. lC.4l for the median 
extinction corrected values, as expected. We therefore conclude 
that contamination driven by foreground extinction cannot ac¬ 
count for the observed trend presented in Fig. [3] 


We note that our requirement of a PACS 70 pm detection 
and the inclusion of longer wavelengths covering the peak of 
the cold-dust SE D, in combination wi th meticulous short wave¬ 
length selection (iMegeath et al.ll20l^ . virtually guarantees the 
presence of an envelope. Therefore we do not consider Class I 
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Fig. C.4. Dependence of the fraction of Class 0 protostars on 
Tboi misclassification due to foreground extinction. Diamonds 
indicate the raw uncorrected observed fractions. Squares indi¬ 
cate the fractions derived by correcting for the sample median 
Ay extinction levels toward protostars in each region. X-symbols 
indicate the fractions derived by correcting for the maximum Ay 
observed toward protostars in each region. 


Fig. C.3. Top; Probability that a source with a given Tboi >70 K 
will have Tboi,ext <70 K, marginalizing over a uniform fore¬ 
ground Ay extinction distribution up to the listed value (shown 
in color). There is no Class 0 contamination from Class I models 
with Tboi > 400 K for the Ay values we consider here. Bottom: 
Probability that a source is correctly classified as a Class 0 
source (has Tboi < 70 K). 


protostar contamination from m ore evolved sources to be a sig¬ 
nificant source or error (see also iHeiderman & Evanil2015h . 
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